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HIGHLIGHTS 


•  A  dynamic  thermal-hydraulic  model  is  developed  for  different  stack  flow  patterns. 

•  The  flow  rate  inhomogeneity  and  reversible  entropic  heat  are  included  in  the  model. 

•  The  model  is  benchmarked  by  experiments  with  mean  relative  error  of  5.85%. 

•  The  battery  performance  is  significantly  influenced  by  the  stack  flow  patterns. 

•  The  serpentine-parallel  pattern  is  optimal  to  enhance  the  battery  performance. 
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The  present  study  focuses  on  dynamic  thermal-hydraulic  modeling  for  the  all-vanadium  flow  battery  and 
investigations  on  the  impact  of  stack  flow  patterns  on  battery  performance.  The  inhomogeneity  of  flow 
rate  distribution  and  reversible  entropic  heat  are  included  in  the  thermal-hydraulic  model.  The  elec¬ 
trolyte  temperature  in  tanks  is  modeled  with  the  finite  element  modeling  (FEM)  technique  considering 
the  possible  non-uniform  distribution  of  electrolyte  temperature.  Results  show  that  the  established 
model  predicts  electrolyte  temperature  accurately  under  various  ambient  temperatures  and  current 
densities.  Significant  temperature  gradients  exist  in  the  battery  system  at  extremely  low  flow  rates,  while 
the  electrolyte  temperature  tends  to  be  the  same  in  different  components  under  relatively  high  flow 
rates.  Three  stack  flow  patterns  including  flow  without  distribution  channels  and  two  cases  of  flow  with 
distribution  channels  are  compared  to  investigate  their  effects  on  battery  performance.  It  is  found  that 
the  flow  rates  are  not  uniformly  distributed  in  cells  especially  when  the  stack  is  not  well  designed,  while 
adding  distribution  channels  alleviates  the  inhomogeneous  phenomenon.  By  comparing  the  three  flow 
patterns,  it  is  found  that  the  serpentine— parallel  pattern  is  preferable  and  effectively  controls  the  uni¬ 
formity  of  flow  rates,  pressure  drop  and  electrolyte  temperature  all  at  expected  levels. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Although  fossil  fuels  are  not  a  clean  energy  resource  because  of 
the  serious  pollutant  emission,  coal-fired  energy  has  played  a 
dominant  role  in  the  electrical  power  generation  industry  for  a  long 
time  [1,2].  Renewable  energy  has  been  viewed  as  the  promising 
substitute  for  coal-fired  energy  during  the  past  years  because  it  can 
save  the  limited  fossil  fuels  while  reducing  pollution  and  green¬ 
house  gas  emissions.  The  intermittent  nature  of  renewable  energy, 
however,  has  posed  a  rigorous  challenge  for  widespread  applica¬ 
tion.  Recent  advances  and  developments  in  large-capacity  batteries 
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have  made  it  possible  for  renewable  energy  to  be  used  in  a 
continuous  and  steady  way  [3,4].  Among  the  multitude  of  large- 
capacity  batteries,  the  all-vanadium  redox  flow  battery  (VRB) 
proposed  by  Skyllas-Kazacos  and  co-workers  [5,6]  in  the  University 
of  New  South  Wales  has  been  viewed  as  the  most  promising.  VRBs 
have  also  been  successfully  commercialized  in  many  countries  for 
the  unique  advantages  of  elimination  of  cross  contamination,  in¬ 
dependent  capacity  and  output  power  design,  tolerance  to  deep 
discharge,  high  energy  efficiency  and  long  life  cycle  [6—9]. 

Thermal  modeling  is  a  valuable  tool  that  should  be  considered 
for  VRB  design  and  operation.  For  one  reason,  the  electrolytes 
should  be  strictly  controlled  within  a  specific  temperature  range  for 
the  battery  to  work  efficiently  10].  For  instance,  for  typical  elec¬ 
trolytes  with  1.8-2  M  vanadium  sulfate  in  2.5—3  M  sulfuric  acid, 
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Nomenclature 

AS® 

standard  molar  reaction  entropy  (J  K-1  mol-1) 

T 

temperature  (I<) 

A 

surface  or  cross-section  area  (m2) 

t 

time  (h) 

a 

sign  factor 

V 

volume  (L) 

cp 

specific  heat  of  electrolyte  (J  kg-1  I<-1) 

w 

width  (m) 

c 

concentration  of  species  (mol  IT1) 

X 

characteristic  length  (m) 

d 

diameter  (m) 

z 

number  of  moles  of  electrons  exchanged  in  a  reaction 

Dh 

hydraulic  diameter  (m) 

df 

mean  fiber  diameter  of  electrode  (m) 

Greek 

symbol 

E® 

standard  electrode  potential  (V) 

P 

dynamic  viscosity  (Pa  s) 

F 

Faraday’s  constant  (C  mol-1) 

e 

porosity  of  graphite  felt  electrode 

f 

friction  resistance  coefficient 

6 

thickness  (m) 

g 

acceleration  due  to  gravity  (m  s-2) 

P 

density  (kg  m-3) 

Gr 

Grashof  number 

V 

volumetric  thermal  expansion  coefficient 

AC® 

standard  molar  Gibbs  free  reaction  enthalpy  (J  mol-1) 

J 

activity  coefficient 

AH® 

standard  molar  reaction  enthalpy  (J  mol-1) 

h 

overall  heat  transfer  coefficient  (W  m-2  K-1) 

Subscript 

L 

length  (m) 

c 

distribution  channels 

I 

current  (A) 

com 

flow  combination 

Kp 

permeability  coefficient 

e 

electrode 

I<ck 

Carman-Kozeny  constant 

ec 

electrode  with  channels  on  surface 

I< 

resistance  coefficient 

f 

fluid 

k 

thermal  conductivity  (W  m-1  K-1) 

i 

inner  surface 

N 

number  of  cells 

in 

inlet 

Nu 

Nusselt  number 

m 

manifold  channel 

P 

power  (W) 

mi 

input  manifold  channel 

P 

pressure  (Pa) 

mo 

output  manifold  channel 

Ap 

pressure  drop  (Pa) 

0 

outer  surface 

Pr 

Prandtl  number 

P 

pipe 

Q 

operating  flow  rate  (cm3  s-1) 

Pi 

pipe  connecting  the  inlet  of  stack  and  outlet  of  tank 

Q 

flow  rate  in  cells  (cm3  s-1) 

po 

pipe  connecting  the  outlet  of  stack  and  inlet  of  tank 

R 

gas  constant  (J  I<-1  mol-1) 

r 

rib 

r 

radius  (m) 

sep 

flow  separation 

Ra 

Rayleigh  number 

s 

stack 

Re 

Reynolds  number 

t 

tank 

Rc 

stack  resistance  during  charge  (Q) 

w 

wall 

Ra 

stack  resistance  during  discharge  (Q) 

+ 

positive  side 

SOC 

state  of  charge  (%) 

— 

negative  side 

thermal  precipitation  of  V5+  happens  when  the  temperature  is 
higher  than  40  °C  over  extended  periods,  while  precipitation  of 
V2+/v3+  occurs  when  the  electrolyte  temperature  drops  below  5  °C 
[11].  Battery  efficiency  cannot  be  confirmed  and  the  precipitation 
can  even  block  the  electrolyte  channels  if  the  temperature  is  left 
uncontrolled.  Although  precipitation  can  be  minimized  by  lowering 
the  vanadium  ion  concentration  below  1.8  M,  this  leads  to  a 
reduced  energy  density  that  is  often  undesirable.  For  this  purpose,  a 
temperature  control  system  is  highly  desirable  for  the  high- 
efficiency  operation  of  the  VRB,  although  such  a  battery  control 
system  can  only  be  established  based  on  the  accurate  on-line  pre¬ 
diction  of  electrolyte  temperature.  Furthermore,  accurate  temper¬ 
ature  prediction  serves  as  the  requisite  for  modeling  the  electrical 
behavior  of  the  VRB,  which  is  vital  for  interfacing  the  battery  with 
other  power  electronic  devices  [12].  This  is  because  of  many  crucial 
parameters  in  electrical  models,  for  example  electrode  potential, 
cell  resistance  and  SOC  determination  are  all  related  to  the  elec¬ 
trolyte  temperature  [13]. 

Though  of  great  importance,  studies  on  thermal  modeling  of 
VRB  have  been  considerably  limited.  A  two-dimensional  numerical 
model  was  established  by  Al-Fetlawi  et  al.  14].  Such  a  numerical 
model  was  able  to  acquire  detailed  temperature  distribution  in¬ 
formation  in  the  battery  and  is  helpful  for  the  battery  design.  The 


high  calculation  complexity  and  long  simulating  times,  however, 
was  not  suitable  for  the  dynamic  control  of  electrolyte  temperature. 
Tang  et  al.  [15  first  proposed  a  dynamic  thermal  model  considering 
the  energy  and  mass  conservation.  The  electrolyte  temperatures 
under  both  constant  and  varying  surrounding  temperature  were 
studied.  Further,  the  effect  of  self-discharge  reactions  was  added  to 
the  model  to  acquire  a  more  reasonable  solution  of  temperature 
[16].  One  limitation  was  that  the  reversible  entropic  heat  [20]  was 
simplified  and  incorporated  into  a  cell  resistance  term.  In  more 
recent  papers,  Tang  et  al.  [17,18]  modeled  the  effect  of  shunt  current 
and  flow  rate  on  battery  efficiency,  stack  temperature,  pressure 
losses,  etc.  In  the  earlier  work,  a  dynamic  model  was  built  by 
considering  the  impacts  of  pump  energy  losses  on  the  electrolyte 
temperature  [19]  but  several  issues  were  still  expected  to  be 
improved.  The  entropic  heat  from  chemical  reaction  was  consid¬ 
ered,  but  this  heat  source  was  simplified  as  enthalpy  change  of  the 
electrode  reaction,  which  reduced  the  accuracy  of  the  simulations. 
Another  essential  limitation  of  the  earlier  dynamic  models  was  that 
the  effect  of  stack  flow  patterns  was  not  considered  and  flow  rate 
was  assumed  to  be  uniformly  distributed  in  each  cell,  which  actu¬ 
ally  deviated  from  the  practical  situation.  Therefore,  the  earlier 
dynamic  models  are  not  adequate  to  be  used  in  the  practical  stack 
design  because  significant  non-uniform  distribution  of  flow  rates 
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can  occur  when  the  stack  is  not  well  designed.  In  addition,  in  the 
earlier  models,  the  temperature  was  also  considered  to  be  uniform 
in  each  component  of  the  system  which  was  not  always  the  case.  All 
of  these  assumptions  will  affect  the  accuracy  of  the  model  and 
therefore  requires  further  improvement. 

As  discussed  above,  stack  flow  pattern  design  is  crucial,  not  only 
because  it  can  influence  the  accuracy  of  the  thermal  model,  but  also 
because  the  stack  flow  pattern  itself  can  have  an  impact  on  pressure 
drop,  cell  resistance  and  overall  battery  performance.  Similar  to  the 
conclusion  that  different  flow  field  designs  influence  the  flow  dis¬ 
tribution  on  a  specific  graphite  plate  [21],  the  stack  flow  pattern 
design  also  determines  the  flow  rate  distribution  among  different 
cells.  Considering  that  uniform  electrolyte  distribution  in  the  whole 
reaction  region  is  the  prerequisite  of  high-efficiency  operation  of 
the  battery,  stack  flow  pattern  is  an  important  parameter  for 
investigation  and  optimization. 

Many  stack  flow  patterns  have  been  proposed  and  tested  to 
improve  the  efficiency  of  VRB  during  the  past  years.  Apart  from  the 
conventional  flow-through  pattern,  other  stack  flow  patterns  like 
the  serpentine  and  parallel  patterns  proposed  for  use  in  fuel  cells, 
have  also  been  adapted  and  reported  for  the  VRB  [22-24].  Xu  et  al. 
[25]  presented  a  three-dimensional  model  to  compare  the  battery 
efficiencies  under  different  stack  flow  patterns.  The  results  showed 
that  the  serpentine  pattern  is  better  than  the  parallel  and  con¬ 
ventional  flow-through  pattern  considering  the  lower  over¬ 
potential.  In  another  study  [21],  the  flow  rate  distribution  in 
different  channels  of  the  same  plate  was  studied  with  different  flow 
patterns  considered.  These  works  succeeded  in  revealing  how  the 
stack  flow  pattern  design  can  influence  battery  performance  but 
were  only  limited  to  lab-scale  single  cells.  The  effect  of  stack  flow 
pattern  on  the  whole  stack  with  several  parallel  aligned  cells  was 
not  investigated.  Furthermore,  the  thermal  behavior  of  the  battery 
under  different  stack  flow  patterns  was  also  ignored. 

In  the  present  study,  a  dynamic  thermal-hydraulic  model  is 
established.  The  inhomogeneity  of  flow  rate  distribution  and 
reversible  entropic  heat  are  included  in  the  model.  Finite  element 
modeling  (FEM)  is  applied  for  the  tank  to  avoid  errors  caused  by 
non-uniform  temperature  distribution  when  low  flow  rate  and 
poor  mixing  is  applied.  Furthermore,  the  flow  features  in  stack 
under  three  different  flow  patterns,  including  flow  pattern  without 
distribution  channels,  serpentine  pattern  and  serpentine-parallel 
pattern,  are  investigated  and  compared.  The  effects  of  different 
flow  patterns  on  flow  rate  distribution,  pressure  drop  and  elec¬ 
trolyte  temperature  are  clarified.  The  present  study  can  provide  a 
feasible  dynamic  model  for  the  battery  temperature  control  system 
and  can  be  further  coupled  into  electrical  models  of  the  VRB.  By 
integrating  the  hydraulic  calculations,  the  model  is  also  instructive 
for  practical  stack  design  to  control  flow  rate  uniformity,  pressure 
drop  and  electrolyte  temperature  at  expected  level.  The  analysis  of 
flow  patterns  also  gives  deep  insights  into  the  high-efficiency 
design  and  operation  of  VRB. 


converted  to  heat  and  added  to  the  electrolyte.  The  following 
equation  can  be  used  to  evaluate  the  pump  power: 

^pump  —  Ap  x  Q.  (1) 

where  Ap  is  the  pressure  drop  in  the  stack  and  pipe  circuit  caused 
by  hydraulic  resistance,  Q.  is  the  operating  flow  rate. 

2.1.  Pressure  drop  in  stack 

The  pressure  drop  in  the  stack  comes  mainly  from  the  hydraulic 
resistances  of  the  flow  fields  in  the  cells  and  manifold  channels  that 
connect  the  cells. 


2.1.1.  Pressure  drop  in  cells 

The  stack  consists  of  several  parallel  aligned  cells  with  internal 
flow  fields.  Different  electrolyte  flow  patterns  in  the  cells  can  result 
in  distinctly  different  flow/ion  distribution  and  pressure  drop.  In  the 
present  study,  all  the  flows  are  assumed  laminar  considering  the 
relatively  low  flow  rate  during  normal  operation  of  the  VRB.  This 
assumption  is  also  confirmed  in  the  calculation.  In  this  paper,  three 
different  stack  flow  patterns  are  considered,  including  a  flow  pattern 
without  distribution  channels  in  the  bipolar  plate,  a  serpentine 
pattern  and  a  serpentine-parallel  pattern,  as  illustrated  in  Fig.  1. 

Case  1 :  Flow  without  distribution  channels 


This  flow  pattern  corresponds  to  the  first  design  in  Fig.  1.  In  this 
flowing  pattern,  the  electrolyte  flows  directly  through  the  porous 
electrode.  According  to  the  Darcy’s  law,  the  pressure  drop  can  be 
described  as: 


Ape 


pieQ 

l(pAe 


(2) 


where  Kp  is  the  permeability  coefficient  that  can  be  calculated  ac¬ 
cording  to  the  Kozeny-Carman  equation  as  below: 


dje 


where  df  is  the  mean  fiber  diameter  of  the  electrode,  e  is  the 
porosity  of  graphite  felt  electrode,  I<Ck  is  the  Carman-Kozeny  con¬ 
stant  for  the  fibrous  medium. 


Case  2:  Flow  with  distribution  channels 

The  two  flow  patterns  correspond  to  the  second  and  third  ones 
in  Fig.  1.  The  electrolytes  flow  through  the  designed  channels  while 
meanwhile  permeating  into  the  porous  electrodes.  Similarly,  the 
pressure  drop  through  the  porous  electrode  can  be  calculated  by: 


2.  Pressure  drop  and  pump  power 

In  the  VRB  system,  the  motor  driven  pump  provides  mechanical 
power  to  overcome  the  pressure  drop  (hydraulic  resistances)  and 
sustain  the  electrolyte  flow  in  the  positive  and  negative  half-cell 
electrolyte  loops.  Part  of  the  total  electrical  energy  input  into  the 
motor-pump  system  converts  to  the  mechanical  power  mentioned 
above.  The  remaining  part  of  the  electrical  energy  input  cannot  be 
used  by  the  system  and  is  converted  to  heat  that  is  dissipated  to  the 
surroundings.  It  should  be  noted  that  the  pump  power  in  the  pre¬ 
sent  study  is  defined  as  the  mechanical  power  portion  of  the  total 
electrical  energy  input.  From  the  view  of  energy  conservation,  the 
continuous  pump  power  supplied  to  the  electrolytes  is  eventually 


Apec 


pLcQe 


KA* 


(4) 


where  Qe  is  the  flow  rate  in  porous  electrodes  of  each  cell,  Lc  is  the 
length  of  distribution  channel  and  Aec  is  the  characteristic  area  of 
electrode  with  distribution  channel  on  the  surface,  which  can  be 
expressed  by: 

Aec  =  (wc  +  wr)6e  (5) 


where  wc  and  wr  are  respectively  the  widths  of  distribution  channel 
and  rib,  6e  the  thickness  of  the  porous  electrode.  The  pressure  drop 
through  the  distribution  channels  has  basically  two  sources.  One  is 
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Flow  pattern  without 
distribution  channels 


Serpentine  pattern 


.Distribution 

channel 


Rib 


Serpentine-parallel 

pattern 


Fig.  1.  Three  stack  flow  patterns  considered  in  the  present  study. 


the  friction  pressure  drop  caused  by  the  shear  stress  from  the  wall- 
fluid  interface  and  within  the  fluid  due  to  the  viscosity  of  the  fluid. 
The  other  source  is  the  form  pressure  drop  caused  by  sudden 
change  of  fluid  direction  or  velocity.  Taking  both  parts  into  account, 
the  pressure  drop  through  the  channel  is  described  as: 


where  fc  is  the  friction  resistance  coefficient  of  the  distribution 
channels  and  can  be  calculated  by: 


f  _  68.36  fiAc  ^ 

fc  pDhcQc  (  ] 

where  the  value  of  68.36  is  decided  by  the  specific  geometry  of 
channels  [26].  Under  steady  flow  state,  A pc  is  equal  to  A pec  and  the 
values  of  Qe  and  Qc  can  be  obtained.  Under  steady  state  flow,  Apc  is 
equal  to  A pec  and  the  values  of  Qe  and  Qc  can  be  obtained.  Assuming 
that  Qc  =  mQ,  the  pressure  drop  in  cells  can  finally  be  expressed  by: 


Inlet  manifold  channels 


Fig.  2.  Hydraulic  network  representing  the  electrolyte  flow  in  stack. 


(8) 


fmLm\ 


P 

2  Al 


(10) 


2.1.2.  Pressure  drop  in  the  manifold  channels 

The  resistance  through  the  inlet  and  outlet  manifold  channels 
between  two  adjacent  cells  is  the  sum  of  resistance  of  rough  pipe 
flow  and  the  resistance  due  to  the  separating  and  combining  of 
flow,  as  expressed  by: 


Rmi  = 


fmLm\ 

) 


P 

2  A2 


(9) 


where  I<sep  and  I<c om  are  respectively  the  resistance  coefficients  of 
separating  and  combining  of  flow  [27],  Lm  is  the  distance  between 
two  adjacent  cells. 


2.1.3.  Pressure  drop  through  the  whole  stack 

Electrolyte  flow  in  the  whole  stack  can  be  treated  as  a  hydraulic 
network  where  parallel  and  tandem  flows  happen  simultaneously 
[28  ,  as  shown  in  Fig.  2.  As  can  be  seen  in  Fig.  2,  when  the  elec¬ 
trolytes  flow  along  the  circuit  represented  by  the  solid  red  line  (in 
the  web  version),  the  pressure  drop  can  be  expressed  as: 


Pijc  Po,k  ~ 


=  K< 


fm^n 


vsep  ' 


Kcom  + 
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fmLm  \  P 


Dhm  )  2A^ 


^  \2  (  ^  \2  fv- 

)  +  (  X  q> )  +  ■■■+(  X'J* 

i  =  k  J  \i  =  k+ 1  /  \i  =  n 

"  \2  /  "  \2  (n  ^ 

X'Zi)  +[  X  qi  )  +-+(X‘7i 

i  =  k  /  \ /  =  /<+!  /  \i  =  n 
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Qn 
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(Without  distribution  channels) 


fmlm\  P 


v  fmLm\  P 

Acom  ~r 


Dhm  7  2A2, 
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X* 

i=k 


X® 
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X* 
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+  (/cDhCc  +  Kform)  2 Af" 


(With  distribution  channels) 


(12) 
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where  q,  represents  the  flow  rate  in  the  i-th  cell.  When  the  elec¬ 
trolytes  flow  along  the  path  represented  by  the  doted  blue  line  (in 
the  web  version),  the  pressure  drop  can  be  expressed  as: 


PiJ<  Po,k 


HLe 

KpAe 


clk  1  ■ 


u__ 

Dhc 


K, 


form 


P  2 

2A&k~v 


(Without  distribution  channels) 


(With  distribution  channels) 


(13) 


By  solving  Eqs.  (11 )— (13),  the  flow  rate  through  the  fc-th  cell  can 
be  obtained  as: 


Qk- 1 


Qn 


pKpAe 
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I<C  C 
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2 fmLn 
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2/lLeA^ 
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(Without  distribution  channels) 
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^ _ Uhm  J 
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(With  distribution  channels) 


(14) 


The  above  equations  can  be  solved  by  applying  an  iterative 
method,  after  which  the  flow  rate  distributions  among  the  cells  can 
be  determined.  Then  the  pressure  drop  through  the  whole  stack 
can  be  calculated  by  substituting  q \  into  Eqs.  (11)  and  (12). 

2.2.  Pressure  drop  in  the  pipe  circuit 


Q+cpP(Tpi  —  Ts)  +  Q-Cpp(Tpi  —  Ts)  +  I2Rd  +  Ps 


(15) 


Apart  from  the  pressure  drop  caused  by  the  flow  in  battery  stack, 
the  pressure  drop  also  occurs  when  the  electrolytes  flow  through  the 
pipe  circuit.  Similar  to  the  flow  in  the  distribution  channels,  the 
pressure  drop  in  the  pipe  circuit  can  also  be  divided  into  the  friction 
pressure  drop  generated  along  straight  pipes  with  rough  surfaces 
and  the  form  loss  at  elbows,  valves,  inlet  and  outlet  of  electrolyte 
reservoirs.  All  these  pressure  drops  can  be  calculated  according  to 
the  discussion  in  Section  2.1  and  will  not  be  further  repeated  here. 

3.  Dynamic  thermal-hydraulic  model 

The  dynamic  thermal-hydraulic  model  is  established  based  on 
the  energy  conservation  inside  the  stack,  tanks  and  pipe  circuits. 
The  following  assumptions  are  made  in  the  present  study: 

(1)  The  electrolyte  flow  is  incompressible; 

(2)  The  positive  and  negative  electrolytes  do  not  change  in  vol¬ 
ume  during  operation; 

(3)  No  side  reactions  occur  when  state  of  charge  is  maintained 
within  the  range  10%— 90%. 


where 

Ts :  the  temperature  of  stack, 

Tpi:  the  temperature  of  pipe  between  stack  inlet  and  tank  outlet, 
Rd :  the  stack  resistance  during  discharge, 

Ed:  the  open  circuit  voltage  for  electrode  reaction  during 
discharge, 

Ps:  the  pump  power  through  the  stack. 

The  first  two  terms  on  the  right  hand  side  of  Eq.  (15)  are  the 
heat  accumulations  due  to  the  electrolytes  exchange  between 
pipe  and  stack.  The  third  term  represents  the  ohmic  resistance 
heat  which  is  positive  during  both  charge  and  discharge.  The 
fourth  term  is  the  pump  power  consumed  in  the  stack,  as  have 
been  discussed  in  Section  2.  The  last  term  represents  the 
reversible  entropic  heat  and  can  be  periodically  negative  and 
positive  during  charge  and  discharge.  As  the  entropic  heat  can 
significantly  influence  the  thermal  behavior  especially  for  the 
temperature  fluctuation  during  charge  and  discharge,  it  is 
included  as  a  separate  term  in  the  present  model.  The  entropic 
heat  can  be  obtained  by  introducing  the  Nernst  Equation  as 
follows: 


3.1.  Thermal  behavior  in  stack 

During  discharge,  the  energy  balance  equation  of  the  stack  is 
expressed  as  below: 


E  =  £e+Ein 
zF 


( cvo+cv2+  ch  A  f  Tvo+TV2+Th+\ 

\  cvo2+cv3+  /  V  ^vo2+^v3+  J 


(16) 
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where  E®  is  the  standard  electrode  potential  for  electrode  reaction 
when  the  activity  of  each  product  and  reactant  is  unity.  The  product 
term  containing  y  can  be  further  set  equal  to  1  because  the  y\  terms 
can  be  assumed  to  approximately  cancel  each  other.  According  to 
the  chemical  thermodynamics,  the  following  equation  is 
established: 

e  AC0  AH0 -IAS0 

E  =-^r  = - zf —  (17) 

where  AG0,  AH0,  and  AS0  are  the  molar  Gibbs  free  reaction 
enthalpy,  molar  reaction  enthalpy  and  molar  reaction  entropy  at 
standard  conditions  respectively  and  the  values  are  listed  in  Table  1 
[29].  The  value  of  dE/dT  can  be  obtained  by  substituting  (17)  into 
(16),  and  the  entropic  heat  can  be  expressed  as: 


IT—  =  — 

dr  zF 


JS0  +  R  In 


:VO+  Cv2+  CH 


(18) 


The  vanadium  concentrations  inside  the  cells  change  with  the 
process  of  electrochemical  reactions,  resulting  in  different  con¬ 
centrations  at  the  inlets  and  outlets  of  the  cells.  Therefore,  the 
concentrations  of  species  in  Eq.  (18)  are  referred  to  the  average 
values  of  inlet  and  outlet.  The  input  concentration  is  equal  to  the 
concentration  in  tank  and  can  be  expressed  as: 


Cin(t)  =  initial  +  ^  |  ^dt  (19) 

0 

where  Vt  is  the  volume  of  tank,  a  is  the  sign  factor  the  value  of 
which  is  1  for  products  while  -1  for  the  reactants.  The  output 
concentration  relates  to  the  flow  rate  applied  to  the  system  and  is 
expressed  as: 


Cout(t)  =  Qn(t)+=^ff  (20) 

Above  is  the  description  of  energy  balance  equation  during 
discharge.  Similarly,  the  energy  conservation  during  charge  is 
described  as  follow: 


CppVs 


dTs 
d  t 


Q+cpp(Tpi  —  Ts)  +  Q.-Cpp(Jpi  —  Ts )  +  I2Rc  +  Ps 


IT  d£c 

lTsd T 


(21) 


where  Rc  represents  the  ohmic  resistance  during  charge.  The 
entropic  heat  can  be  deduced  with  the  same  method  as  that  of  the 
discharging  process. 


3.2.  Thermal  behavior  in  tanks 


The  temperature  of  the  electrolytes  in  the  tanks  is  determined 
by  the  heat  transfer  between  the  pipe  and  tank  brought  by  flow  and 
the  heat  dissipated  to  the  surroundings.  In  the  present  study,  FEM  is 
implemented  to  tackle  the  possible  non-uniformed  distribution 
phenomenon  where  imperfect  mixing  occurs  in  the  tanks. 
Considering  the  symmetry  of  both  tank  geometry  and  flow,  one¬ 
dimensional  FEM  along  the  axial  direction  of  tank  is  conducted  to 
improve  modeling  accuracy  when  meanwhile  maintain  instanta- 
neity.  The  energy  balance  equation  reads: 


Table  1 

Thermodynamic  data  for  species  in  redox  reaction  at  temperature  of  298.15  K. 


Formula 

State 

AG®(kJ  mol-1) 

AH®(kJ  mol1) 

AS®(J  mol  1  K-1) 

V2+ 

aq 

-218.0 

-226.0 

-130.0 

V3+ 

aq 

-251.3 

-259.0 

-230.0 

V02+ 

aq 

-446.4 

-486.6 

-133.9 

VOJ 

aq 

-587.0 

-649.8 

-42.3 

H4* 

aq 

0 

0 

0 

h2o 

aq 

-237.2 

-285.8 

69.9 

r  dT 

cpPVt, l  ^  =  Q+CpP(7po  —  Tt) i)  +  7it  iv4t  i  (Tair  —  Tt  i) 

CpPVt, =  Q.+cpP(7t,l  -  Tt, 2)  +  ^t,2^t,2  (Tair  -  Tt, 2) 

dr  •  (22) 

cpPvtj~ =  Q+cPp(rtj-i  -  Tfj)  +  htJAtj(Tair  -  Ttj ) 

CppVt  i-^-  =  Q+Cpp(Tt  i_ 1  -  Tt  i)  +  htjAt  i(Tair  -  Ttj ) 

where  Ttj  is  the  temperature  of  electrolytes  in  the  j-th  volume 
element,  Vtj  is  the  volume  of  the  j-th  volume  element,  htj  and  Atj 
are  respectively  the  overall  heat  transfer  coefficient  and  superficial 
area  of  tank  that  encompass  the  j-th  volume  element. 

3.3.  Thermal  behavior  in  pipes 

Similar  to  the  energy  balance  of  tanks,  the  temperatures  of 
electrolytes  in  pipes  are  controlled  by  the  following  equations: 

CppVp0  it  Q+Cpp(Ts  —  Tp0 )  +  hpoApo  (Tair  —  Tpo)  (23) 

dT 

CpPVpi  ^  =  Q+Cpp(Tt  i  —  Tpj)  +  hpiApi  (Tair  —  Tpj)  (24) 

where  Tpo  denotes  the  temperature  of  pipe  that  connect  outlet  of 
stack  and  inlet  of  tank. 

3.4.  Calculation  of  heat  dissipation 

The  outer  surfaces  of  the  tanks  and  pipes  are  supposed  to  be 
exposed  to  the  natural  convection  environment.  Also,  the  inner 
surfaces  of  tanks  can  also  be  treated  as  natural  convection  pro¬ 
cesses  because  the  flow  velocity  in  the  tanks  can  be  approximated 
to  be  low.  Instead,  forced  convection  heat  transfer  happens  be¬ 
tween  the  inner  walls  of  the  pipes  and  electrolytes  for  higher  flow 
velocities.  The  convective  heat  transfer  coefficient  can  be  calculated 
by: 


Table  2 

Empirical  coefficients  of  experimental  correlation  for  different  geometry. 


Geometry 

B 

c 

Condition 

Vertical  plate  or  cylinder 

0.59 

1/4 

104  <  Gr  <3  x  109 

Horizontal  plate  (high  temperature  upward) 

0.54 

1/4 

104  <Ra<  107 

Horizontal  plate  (high  temperature  upward) 

0.15 

1/4 

107  <Ra<  1011 

Horizontal  plate  (high  temperature 
downward) 

0.27 

1/4 

105  <Ra<  1010 

Z.  Wei  et  al.  /  Journal  of  Power  Sources  260  (2014)  89-99 


95 


h 


conv 


(25) 


where  Nu  is  the  Nusselt  number,  kf  the  thermal  conductivity  of 
fluid,  x  the  characteristic  length.  The  Nusselt  number  for  the  natural 
convection  heat  transfer  problem  can  be  acquired  by  the  following 
experimental  correlation: 

Nu  =  b(Ra)c  (26) 


where  b  and  c  are  the  empirical  coefficients  which  are  listed  in 
Table  2  for  different  geometry  [26  ,  Ra  is  the  Rayleigh  number  and 
can  be  expressed  by  the  product  of  Grashof  number  (Gr)  and 
Prandtl  number  (Pr).  These  two  numbers  can  be  calculated  by: 


Pr  = 


cjm 


(27) 


Gr 


gpjffhTx3 


(28) 


where  cp/is  the  specific  heat  of  the  fluid,  /if  the  dynamic  viscosity  of 
fluid,  1 6  the  volumetric  thermal  expansion  coefficient,  AT  the  tem¬ 
perature  difference  between  wall  and  fluid. 

For  forced  convection  heat  transfer  inside  the  pipes,  the  Nusselt 
number  can  be  calculated  by  the  Sieder-Tate  equation  [26]: 


Nu  =  1.86 


/ RePr\ i 
\LP/  dp) 


(29) 


where  Re  represents  the  Reynolds  number,  r\/r\w  describes  the  in¬ 
fluence  of  temperature  gradient  in  pipes  and  is  empirically  decided 
as  0.95  in  engineering  problems. 

After  the  convective  heat  transfer  coefficients  are  obtained,  the 
overall  heat  transfer  coefficients  can  be  determined  by: 

(  1 

7 - - - j - - - ,  (Cylindrical  wall) 

l  r  ,  r  +  C7W  r 

7 - — b  i —  In - b  t - — z~r 

.  hi  kw  r  (r  +  0w)ho 

i 

1 — K — V  (Flat  wall) 

hi  kw  h0 

(30) 

where  hi  and  h0  denote  the  convective  heat  transfer  coefficients  of 
the  inner  and  outer  surfaces,  respectively. 

4.  Results  and  discussion 

4.1.  Simulation  parameters 


Table  3 

Parameters  adopted  in  simulation. 


Parameters  Value 


Number  of  cells 
Dimensions  of  stack 
Area  of  electrode 

Thickness  of  graphite  felt  electrode 
Porosity  of  felt  electrode 
Fiber  diameter 
Kozeny— Carman  constant 
Width  of  distribution  channels 
Depth  of  distribution  channels 
Diameter  of  manifold  channels 
Concentrations  of  electrolytes 


Density  of  electrolytes 
Specific  heat  of  electrolytes 
Dynamic  viscosity  of  electrolytes 
Diameter  of  pipes 
Length  of  pipes 

Volume  of  electrolytes  in  each  tank 
Dimensions  of  tank 
Thickness  of  tank  wall 
Thermal  conductivity  of  tank  wall 
Overall  heat  transfer  coefficient 
(top  wall  of  tank) 

Overall  heat  transfer  coefficient 
(bottom  wall  of  tank) 

Overall  heat  transfer  coefficient 
(side  wall  of  tank) 

Overall  heat  transfer 
coefficient  (pipes) 

Stoichiometric  flow  rate  for 

current  density  of  40/100  mA  cm-2 
and  SOC  range  of  10%-90% 

Applied  flow  rate 


14 

0.440  m  x  0.34  m  x  0.2  m 
0.0875  m2 
4  x  10-3  m 
0.68 

3  x  10-15  m 
5.0 

6  x  lO"3  m 
2  x  10~3  m 
12  x  1CT3  m 

Positive:  1.5  M  VOS04+3  M  H2S04 

Negative: 0.75  M  V2(S04)3 

+  2.25  M  H2S04 

1400  kg  nrr3 

3200  J  kg"1  K1 

4.928  x  10  3  Pa  s 

1  x  10'2  m 

3.56  m 

7.4  L 

0.2  m  x  0.2  m  x  0.4  m 
1  x  10-2  m 
0.16  Wm1  K_1 
2.995  W  m  2  K1 

1.703  W  m  2  K_1 

2.873  W  m  2  K1 

3.806  W  m  2  K1 

33.9/84.6  cm3  s"1 


125  cm3  s- 1 


discharging  times  are  determined  by  the  below  Faraday’s  law  of 
electrolysis: 


tc  = 


czFlVt+yf 


Nlr 


Vsoc 


(31) 


fd 


czF(vt  +  f) 

NId 


Vsoc 


(32) 


where  psoc  represents  the  percent  range  of  SOC,  which  is  80%  for 
the  present  simulation.  Although  Tang  et  al.  show  that  ion  diffusion 
across  the  membrane  can  lead  to  considerable  heat  generation 
from  the  self-discharge  reactions,  the  effect  on  electrolyte  tem¬ 
perature  are  shown  to  be  most  significant  when  the  pumps  were 
turned  off  for  several  hours.  In  this  study,  the  pumps  are  assumed  to 
remain  on  during  the  continuous  charge-discharge  cycling,  so  the 
thermal  effect  of  the  self-discharge  reactions  can  be  assumed  to  be 
minimal. 


A  kilowatt  class  VRB  system  with  14  cells  in  the  stack  is  simu¬ 
lated  with  the  model  developed  above.  Parameters  adopted  in  the 
simulation  are  summarized  in  Table  3.  The  electrolytes  flow  directly 
through  the  porous  electrode  without  additional  distribution 
channels.  The  SOC  of  the  battery  is  controlled  between  10%  and 
90%.  The  modeling  is  conducted  for  25  charging— discharging  cy¬ 
cles.  The  dynamic  model,  characterized  by  a  set  of  differential 
equations  which  can  be  treated  as  a  nonlinear  system,  is  solved  in 
the  Matlab/Simulink  environment. 

The  model  assumes  that  no  side  reactions  or  ion  diffusion  across 
the  membrane  occurs  during  operation  and  a  coulomb  efficiency  of 
100%  was  applied.  Under  this  condition,  the  charging  and 


4.2.  Validation  of  model 

To  verify  the  dynamic  thermal-dynamic  model  established  in 
the  present  study,  the  modeling  results  are  compared  with  the 
experimental  data  from  another  published  paper  [30].  Experiments 
are  conducted  under  different  ambient  temperatures.  The  charging 
current  density  is  maintained  as  40  mA  cm-2  while  the  discharging 
current  density  changes  from  30  mA  cnrT2  to  100  mA  cm-2.  The 
flow  pattern  without  distribution  channels  was  applied  to  the  stack 
design.  The  comparison  of  modeled  and  measured  electrolyte 
temperatures  in  stack  under  different  ambient  temperatures  is 
shown  in  Fig.  3.  As  can  be  seen,  the  modeled  electrolyte 
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temperatures  correspond  well  with  the  experimental  ones.  The 
mean  relative  error  (MRE)  and  root  mean  square  error  (RMSE)  of 
modeled  results  are  5.85%  and  0.86  °C,  respectively.  It  can  be 
concluded  that  the  present  dynamic-hydraulic  thermal  model  can 
accurately  predict  the  temperature  of  electrolytes  under  different 
ambient  temperatures  and  different  charging/discharging  current 
densities. 


4.3.  Temperature  responses  for  different  system  components 

Under  a  constant  ambient  temperature  of  25  °C,  the  tempera¬ 
ture  responses  of  different  system  components  for  the  flow  pattern 
without  distribution  channels  are  investigated  in  this  section  using 
the  developed  models.  The  results  are  shown  in  Fig.  4.  For  the  stack 
size  and  electrolyte  volume  used  in  this  simulation,  the  theoretical 
charging  and  discharging  times  are  0.73  and  0.29  h  respectively  at 
the  current  densities  of  40/100  mA  cm-2.  As  can  be  seen,  the 
temperatures  of  all  the  components  increase  gradually  during  the 
first  several  cycles,  then  reach  to  the  steady  state  and  oscillate 
around  the  respective  steady-state  temperature.  Here  the  steady- 
state  temperature  is  defined  as  the  arithmetic  mean  temperature 
of  the  last  charging/discharging  cycle.  It  can  also  be  observed  that 
the  temperatures  of  the  system  components  all  decrease  during 
charge  while  increase  during  discharge  after  the  steady  states  are 
reached.  This  phenomenon  can  be  explained  by  two  reasons.  First, 
the  reaction  rate  of  V5+  reduction  is  kinetically  faster  than  that  of 
V4+  oxidation.  As  a  result,  the  cell  resistance  during  charge  is 
smaller  than  that  during  discharge,  leading  to  more  heat  generated 
during  discharge.  Second,  as  shown  in  Fig.  5,  the  reversible  entropic 
heat  is  periodically  negative  during  charge  and  positive  during 
discharge,  indicating  that  heat  absorption  happens  during  charge 
while  heat  generation  happens  during  discharge. 

The  main  heat  transfer  among  system  components  is  realized  by 
the  flow  of  the  electrolytes,  so  different  flow  rates  can  result  in 
different  heat  transfer  conditions  between  system  components.  As 
can  be  seen  in  Fig.  4,  the  peak  temperature  of  electrolytes  in  the 
stack  is  about  1.6  °C  higher  than  that  in  the  tanks  at  the  flow  rate  of 
5  cm3  s-1.  This  is  due  to  the  fact  that  low  flow  rate  cannot  efficiently 
transfer  heat  from  stack  to  tank,  making  it  possible  for  the  heat 
accumulated  in  the  stack.  Therefore,  long  time  operation  under  low 
flow  rate  should  be  avoided  to  prevent  unexpected  high  tempera¬ 
tures  in  stack,  particularly  when  other  operating  parameters,  for 
example  high  charging/discharging  current  densities  and  ambient 
temperatures,  have  already  increased  the  system  temperature.  By 


Modelled:  40/30  mA  cm'2 
Modelled:  40/100  mA  cm'2 
Experimental:  40/30  mA  cm'2 
Experimental:  40/100  mA  cm"2 


Fig.  3.  Modeled  and  measured  electrolyte  temperature  in  stack  under  different 
ambient  temperatures  (unit:  °C). 


Fig.  4.  Simulation  results  for  the  system  temperature  under  two  different  flow  rates 
(unit:  °C). 


comparison,  the  temperatures  of  system  components  are  approx¬ 
imately  the  same  at  a  flow  rate  of  60  cm3  s-1,  because  the  relatively 
high  flow  rate  effectively  transfers  heat  out  of  the  stack  and  allows 
the  temperature  to  distribute  uniformly  in  the  whole  system. 

4.4.  Influence  of  stack  flow  pattern  design 

In  this  section,  apart  from  the  original  flow  pattern  without 
distribution  channels,  two  other  flow  patterns  with  distribution 
channels  including  the  serpentine  pattern  and  serpentine-parallel 
pattern  are  considered  and  the  effects  of  them  on  battery  perfor¬ 
mance  are  analyzed. 

4.4.1.  Effect  of  manifold  channel  diameter  on  flow  rate  distribution 
among  cells 

The  flow  resistance  of  the  manifold  channel  will  cause  the  non- 
uniform  flow  distribution  among  the  cells.  To  investigate  the  effect 
of  the  manifold  channel  design  and  different  flow  patterns  on  the 
flow  distribution,  a  sensitivity  study  is  carried  out  by  varying  the 
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Fig.  5.  Reversible  entropic  heat  during  25  charging/discharging  cycles. 
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manifold  channel  diameters  for  different  flow  patterns.  The  in¬ 
homogeneity  coefficient  (IC),  defined  as  in  Eq.  (33),  is  adopted  in 
order  to  quantitatively  evaluate  the  inhomogeneity  of  flow  rate 
distribution. 

IC  =  <Jmax  (33) 

Qmin 

where  qmax  denotes  the  maximum  cell  flow  rate  and  qm[n  denotes 
the  minimum  cell  flow  rate.  The  flow  rate  of  125  cm3  s-1  is  adopted 
in  the  sensitivity  study.  The  result  is  shown  in  Fig.  6. 

As  can  be  seen,  the  IC  decreases  with  the  increase  in  manifold 
channel  diameter  for  all  three  patterns.  With  the  same  manifold 
channel  diameter,  the  inhomogeneity  level  for  the  pattern  without 
distribution  channels  is  much  higher  than  the  two  cases  with  dis¬ 
tribution  channels.  Furthermore,  the  IC  of  pattern  without  distri¬ 
bution  channels  is  much  more  sensitive  to  the  manifold  diameter. 
This  can  be  explained  by  analyzing  the  hydraulic  resistance  in  the 
stack.  In  the  case  without  distribution  channels,  the  relatively  low 
flow  resistance  through  the  porous  electrode  is  comparable  to  that 
of  the  manifold  channels.  Therefore,  the  change  in  manifold  resis¬ 
tance  caused  by  changing  diameters  will  significantly  influence  the 
IC  value.  Instead,  for  the  cases  with  distribution  channels,  the 
resistance  in  cell  is  significantly  increased  for  the  existence  of  long 
and  narrow  distribution  channels.  Compared  with  the  cells,  the 
resistance  of  manifold  channels  is  so  small  that  the  change  of  it  has 
no  significant  impact  on  the  IC  value.  Therefore,  the  inhomogeneity 
of  flow  rate  distribution  is  not  much  sensitive  to  manifold  channel 
design. 

From  the  above  analysis,  it  can  be  seen  that  the  non-uniform 
distribution  of  flow  rate  is  significant  if  the  stack  is  not  well 
designed.  Even  if  a  large  manifold  channel  diameter  is  applied,  the 
inhomogeneous  phenomenon  still  exists  for  all  three  flow  patterns. 
Therefore,  the  inhomogeneous  phenomenon  should  be  considered 
in  both  thermal  modeling  and  battery  design.  As  shown  in  Table  3,  a 
manifold  channel  diameter  of  0.012  m  is  implemented  in  the 
following  sections  of  this  paper. 

4.4.2.  Pressure  drop  in  stack 

The  pressure  drop  through  the  stack  is  an  important  parameter 
to  be  observed  because  it  decides  the  pump  power  required  to  flow 
the  electrolytes.  The  pressure  drops  in  the  stack  with  different  flow 
patterns  are  shown  in  Fig.  7.  As  can  be  seen,  the  pressure  drops  with 


Fig.  6.  Inhomogeneity  coefficient  of  flow  rate  with  different  diameters  of  manifold 
channel. 


Fig.  7.  Pressure  drop  in  stack  with  the  three  flow  patterns  under  different  operating 
flow  rates. 


the  three  patterns  all  increase  with  increasing  flow  rate.  The 
serpentine  pattern  yields  the  highest  pressure  drop  under  all  the 
considered  flow  rates,  which  is  due  to  the  high  resistance  of  the 
long  and  narrow  serpentine  channels  used  in  the  simulations.  As 
compared,  the  lower  resistance  of  the  other  two  flow  patterns  helps 
to  keep  the  pressure  drop  at  a  relatively  low  level.  Though  the 
pattern  without  distribution  channels  achieves  lower  pressure  drop 
than  that  of  serpentine-parallel  pattern,  the  difference  is  not 
dramatic. 

It  should  be  mentioned  that  the  present  results  are  not  totally 
consistent  with  the  conclusions  in  other  publication  [25],  where 
the  author  concludes  that  the  pressure  drop  of  flowing  through 
porous  electrode  without  distribution  channels  can  be  higher 
compared  with  serpentine  design  under  low  flow  rate.  This  is 
because  reference  [25]  concentrates  on  a  lab-scale  single  cell  with  a 
small  electrode  area  of  10  cm  x  10  cm,  so  that  the  applied  current 
can  be  considerably  small  and  thus  the  battery  can  be  operated 
under  extremely  low  flow  rate  as  5  mL  s_1.  In  contrast,  the  battery 
in  present  study  is  kilowatt  class.  Though  sometimes  extremely  low 
flow  rate  cannot  be  avoided,  a  relatively  high  flow  rate  is  recom¬ 
mended,  so  a  starting  flow  rate  of  60  cm3  s  1  is  applied  in  this 
analysis.  For  comparison,  the  stoichiometric  flow  rate  calculated 
from  Faraday’s  Law  is  84.6  cm3  s-1  for  a  current  density  of 
100  mA  cm-2  at  an  SOC  of  10%.  That  is  the  reason  that  the  pressure 
drop  for  the  flow  pattern  without  distribution  channels  is  always 
lower  than  that  with  a  serpentine  channel  as  simulated  in  this 
paper. 

It  is  meaningful  to  comprehensively  consider  the  in¬ 
homogeneity  of  flow  distribution  and  pressure  drop  under  the 
three  patterns.  Though  the  serpentine  pattern  facilitates  uniform 
flow  distribution  within  the  cells  under  certain  cell  configurations, 
it  may  also  yield  the  highest  pressure  drop.  This  may  not  be  a  good 
choice  because  a  large  pump  power  is  demanded.  On  the  other 
hand,  for  the  flow  pattern  without  distribution  channels  used  in 
this  study,  the  pressure  drop  is  lower  but  the  flow  rate  distribution 
is  the  most  inhomogeneous,  which  is  also  undesirable.  Only  the 
serpentine-parallel  pattern  produces  a  low  pressure  drop  while 
maintains  uniform  flow  rate  distribution  across  the  cells.  It  should 
be  mentioned  that  two  parallel  channels  are  adopted  in  the  current 
study.  Although  the  pressure  drop  can  be  further  decreased  by 
applying  more  parallel  channels,  this  would  be  at  the  expense  of 
increasing  the  inhomogeneity  of  the  flow  distribution.  Therefore, 
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though  the  serpentine— parallel  pattern  enhances  the  battery  per¬ 
formance  as  a  whole,  designers  still  have  to  make  compromises 
between  the  uniformity  of  flow  rates  and  pressure  drop  to  deter¬ 
mine  the  number  of  parallel  channels  according  to  the  practical 
requirement. 


4.4.3.  Temperature  of  electrolytes 

In  this  section,  the  temperature  of  the  electrolytes  is  discussed. 
As  a  high  flow  rate  of  125  cm3  s-1  is  applied,  the  temperature  in  the 
whole  system  is  uniform,  as  has  been  confirmed  in  Section  4.3. 
Under  a  constant  ambient  temperature  of  25  °C  and  charging/dis¬ 
charging  current  density  of  40/100  mA  cnrT2,  the  temperature  of 
the  electrolytes  was  modeled  and  the  results  under  three  different 
flow  patterns  are  shown  in  Fig.  8.  As  can  be  seen,  the  electrolyte 
temperature  changes  by  following  the  same  trends  although  they 
are  designed  with  different  flow  patterns.  The  steady-state  tem¬ 
peratures  are  37.3  °C,  40.8  °C  and  37.9  °C  respectively  for  the  flow 
pattern  without  distribution  channels,  serpentine  pattern  and 
serpentine-parallel  pattern,  indicating  that  applying  distribution 
channels  can  increase  the  electrolyte  temperature  in  the  current 
case  where  an  increased  pressure  drop  occurs.  It  can  also  be  seen  in 
Fig.  8  that  the  serpentine  pattern  has  contributed  to  increasing  the 
electrolyte  temperature  to  the  thermal  precipitation  region,  which 
could  deteriorate  the  battery  performance  during  operation. 

To  generalize,  more  simulations  were  conducted  under  different 
operating  flow  rates.  Taking  the  case  without  distribution  channels 
as  reference,  the  temperature  increases  caused  by  adding  distri¬ 
bution  channels  under  different  flow  rates  are  shown  in  Fig.  9.  As 
can  be  seen,  the  serpentine  pattern  can  significantly  increase  the 
temperature  of  electrolytes  compared  with  the  case  without  dis¬ 
tribution  channels,  especially  when  a  relatively  high  flow  rate  is 
applied.  For  instance,  a  temperature  increase  of  more  than  5  °C  is 
observed  at  a  flow  rate  of  160  cm3  s-1.  The  significant  temperature 
increase  is  caused  by  two  reasons.  First,  the  high  pressure  drop 
brought  by  the  present  serpentine  pattern  increases  the  pump 
power  consumed  in  the  stack,  which  generates  and  dissipates  more 
heat  to  the  electrolytes.  Second,  the  thermal  dissipation  capacity  of 
the  tank  is  limited  because  a  small  tank  was  applied  in  the  current 
case.  With  a  larger  tank,  the  phenomenon  of  temperature  increase 
can  be  alleviated  to  some  extent  although  it  still  exists.  In  contrast, 
the  serpentine-parallel  pattern  also  leads  to  temperature  increase 
of  electrolytes,  but  the  increase  is  small  enough  to  be  ignored 
especially  at  a  low  operating  flow  rate.  Actually,  the  temperature 
increase  is  also  limited  to  within  1  °C  even  if  the  battery  is  operated 
with  a  high  flow  rate.  This  phenomenon  can  be  explained  by  the 
relatively  low  pressure  drop  increase  caused  by  the  serpentine- 


Fig.  8.  Simulation  result  of  electrolyte  temperatures  under  three  flow  patterns 
(charge/discharge:  40/100  mA  cm-2,  unit:  °C). 


Fig.  9.  Electrolyte  temperature  increase  caused  by  applying  distribution  channels 
(unit:  °C). 


parallel  pattern  used  in  the  present  study  which  can  be  seen  in 

Fig.  7. 

From  the  above  analysis,  it  can  be  seen  that  the  electrolyte 
temperature  of  VRB  is  influenced  by  the  flow  pattern  design. 
Though  the  serpentine  pattern  facilitates  the  uniform  distribution 
of  flow  rates  in  different  cells,  unless  the  channel  dimensions  are 
carefully  selected,  an  increased  pressure  drop  may  arise  leading  to 
increased  electrolyte  temperature  that  can  potentially  cause  the 
thermal  precipitation  of  V5+.  By  comparison,  the  serpentine-par¬ 
allel  pattern  adopted  in  the  present  study  can  significantly  improve 
the  uniformity  of  the  flow  rate  distribution  while  keeping  the 
pressure  drop  and  electrolyte  temperature  close  to  that  of  the 
pattern  without  distribution  channels.  Therefore,  the  choice  of  flow 
pattern  should  be  considered  seriously  during  VRB  design  to 
minimize  the  possibility  of  thermal  precipitation. 

In  Section  4.4,  the  effect  of  flow  pattern  of  stack  on  the  thermal- 
hydraulic  performance  of  battery  is  analyzed  and  the  serpentine- 
parallel  pattern  is  recommended  by  considering  different  factors  to 
enhance  the  battery  performance.  In  practical  applications,  the 
dimensions  of  channels  can  be  further  optimized  for  certain  pur¬ 
pose.  For  instance,  the  pressure  drop  for  the  pattern  without  dis¬ 
tribution  channels  is  sensitive  to  felt  compression  of  the  electrodes. 
Considering  a  battery  with  highly  compressed  felt  electrodes  and 
thus  a  considerably  high  pressure  drop  across  stack,  if  low  pressure 
drop  is  the  major  concern,  inserting  distribution  channels  with 
large  cross-section  areas  can  provide  an  alternative  path  for  the 
electrolyte  flow  and  thus  reduce  the  pressure  drop  compared  with 
the  scenario  without  distribution  channels.  As  a  result,  the  elec¬ 
trolyte  temperature  can  also  be  reduced,  but  a  higher  in¬ 
homogeneity  of  flow  distribution  will  be  the  byproduct.  Therefore, 
the  flow  pattern  and  dimensions  should  be  carefully  evaluated  for  a 
specific  VRB  design  in  the  practical  applications. 


5.  Conclusions 

In  this  study,  the  inhomogeneity  of  flow  rate  distribution  in 
different  cells  is  investigated  and  coupled  in  the  dynamic  thermal- 
hydraulic  model.  Finite  element  modeling  is  implemented  for  the 
tank  to  rule  out  the  possible  error  brought  by  temperature  non¬ 
uniformity.  The  reversible  entropic  heat  is  also  considered  as  a 
heat  source  in  the  stack.  After  the  dynamic  model  was  established, 
three  flow  patterns  are  studied  to  investigate  their  effects  on 
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battery  performance.  The  major  conclusions  are  summarized  as 
follows: 

(1)  The  established  dynamic  thermal-hydraulic  model  accu¬ 
rately  predicts  the  electrolyte  temperature  under  various 
ambient  temperatures  and  charging/discharging  current 
densities.  The  MRE  and  RMSE  of  modeled  results  are  5.85% 
and  0.86  °C,  respectively. 

(2)  The  temperature  distributions  in  the  system  are  distin¬ 
guished  at  different  flow  rate  conditions.  Significant  tem¬ 
perature  gradients  exist  in  the  system  at  an  extremely  low 
flow  rate.  At  a  relatively  high  flow  rate,  instead,  the  electro¬ 
lyte  temperature  distributes  uniformly  in  the  whole  system. 

(3)  Non-uniform  flow  distribution  across  the  cells  in  the  stack 
always  exists  even  when  a  large  manifold  channel  diameter 
is  applied,  while  the  inhomogeneous  phenomenon  can  be 
alleviated  by  adding  distribution  channels  in  the  bipolar 
plate.  The  non-uniform  flow  distribution  phenomenon 
should  therefore  be  considered  in  thermal  modeling  and 
stack  design  of  VRB. 

(4)  The  battery  performance  is  significantly  influenced  by  the 
design  of  stack  flow  pattern.  The  serpentine  pattern  provides 
the  most  uniform  distribution  of  flow  rates,  but  meanwhile 
significantly  increases  the  pressure  drop  and  the  risks  in 
causing  thermal  precipitation  of  V5+  with  the  increased 
electrolyte  temperature.  By  comparison,  the  serpentine- 
parallel  pattern  used  in  the  present  study  combines  the  ad¬ 
vantages  of  the  other  two  flow  patterns  and  effectively 
controls  the  uniformity  of  flow  rate  distribution,  the  pressure 
drop  and  the  electrolyte  temperature. 
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